{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018  by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving the isotropic elastic 2D SH problem by the staggered grid FD method\n",
    "\n",
    "Modelling wave propagation in the solid earth by the acoustic wave equation is a very rough approximation. It is not able to describe S-body waves or Rayleigh/Love surface waves and therefore also not any PS or SP conversions. \n",
    "\n",
    "The most simple elastic approximation is the 2D SH problem, which can be described by partial differential equations mathematically very similar to the 2D acoustic problem. Let's begin with a short review of the 2D SH problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equations of motion for 2D SH wave propagation in an isotropic elastic medium"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described in [this lesson](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/01_Analytical_solutions/4_2D_1D_elastic_acoustic_approx.ipynb), we can decompose the 3D isotropic elastic problem into two decoupled problems, the 2D PSV and 2D SH problem. Assuming only non-zero particle displacements in the y-direction perpendicular to the x-z-plane leads to the **2D SH problem**, where x denotes the horizontal distance and z the depth, wave propagation can be described by the following system of partial differential equations:\n",
    "\n",
    "\\begin{align}\n",
    "\\rho \\frac{\\partial^2 u_y}{\\partial t^2} &= \\frac{\\partial \\sigma_{yx}}{\\partial x} + \\frac{\\partial \\sigma_{yz}}{\\partial z} + f_y, \\notag\\\\\n",
    "\\sigma_{yx} &= \\mu \\frac{\\partial u_{y}}{\\partial x},\\\\\n",
    "\\sigma_{yz} &= \\mu \\frac{\\partial u_{y}}{\\partial z}, \\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "where $\\rho$ is the density and $\\mu$ the shear modulus, $\\rm{u_y}$ the particle displacement in y-direction, $\\sigma_{yx}$, $\\rm{\\sigma_{yz}}$ shear stress components, $f_y$ the directed body force in y-direction, respectively. \n",
    "\n",
    "We can modify this problem by introducing the **particle velocity** $v_y = \\frac{\\partial u_y}{\\partial t}$. Taking the first temporal derivative of the stress-strain relation in eqs.(1) and assuming that the shear modulus $\\mu$ is not time-dependent, yields a system of coupled first order partial differential equations:\n",
    "\n",
    "\\begin{align}\n",
    "\\rho \\frac{\\partial v_y}{\\partial t} &= \\frac{\\partial \\sigma_{yx}}{\\partial x} + \\frac{\\partial \\sigma_{yz}}{\\partial z} + f_y, \\notag\\\\\n",
    "\\frac{\\partial \\sigma_{yx}}{\\partial t} &= \\mu \\frac{\\partial v_{y}}{\\partial x},\\\\\n",
    "\\frac{\\partial \\sigma_{yz}}{\\partial t} &= \\mu \\frac{\\partial v_{y}}{\\partial z}. \\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "This is known as the **stress-velocity formulation** of the 2D SH problem. Next, we solve eqs. (2) by a staggered grid finite difference approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the 2D SH problem by the Finite difference method on a staggered grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As for the previous FD approximation of the 1D and 2D acoustic wave equation, we discretize the elastic equations of motion in time and space on a Cartesian grid. However, compared to the acoustic wave equation, the 2D SH problem eqs. (2) involves three wavefields $v_y,\\; \\sigma_{yx},\\; \\sigma_{yz}$ and two material parameter classes $\\rho,\\; \\mu$. \n",
    "\n",
    "In order to evaluate the spatial derivatives correctly in relation to the wavefields and material parameters, [Jean Virieux](https://jean-virieux.obs.ujf-grenoble.fr/) introduced in [this paper](https://www.isterre.fr/IMG/pdf/gpy_1984_virieux.pdf) the following **Staggered Finite Difference grid** for the 2D SH problem:\n",
    "\n",
    "<img src=\"images/SG_SH-Cart.png\" width=\"50%\">\n",
    "\n",
    "The particle velocity $v_y$ and density $\\rho$ are located at the grid point (i,j), while the shear stress $\\sigma_{yx}$ is displaced by half a grid point in x-direction (i+1/2,j) and shear stress $\\sigma_{yz}$ is shifted by half a grid point in z-direction (i,j+1/2). Furthermore, the shear moduli have to be averaged harmonically:\n",
    "\n",
    "\\begin{align}\n",
    "\\mu_{x} &= \\mu(i+1/2,j) = 2 \\biggl[\\biggl(\\mu^{-1}_{i,j}+\\mu^{-1}_{i+1,j}\\biggr)\\biggr]^{-1} \\notag\\\\\n",
    "\\mu_{z} &= \\mu(i,j+1/2) = 2 \\biggl[\\biggl(\\mu^{-1}_{i,j}+\\mu^{-1}_{i,j+1}\\biggr)\\biggr]^{-1} \\notag\n",
    "\\end{align}\n",
    "\n",
    "$dx,\\; dz$ denotes the spatial distance between two adjacent spatial grid points in x- and z-direction, respectively. The next question is how to approximate the spatial partial derivatives on this staggered grid. \n",
    "\n",
    "We begin with the momentum equation\n",
    "\n",
    "\\begin{equation}\n",
    "\\rho \\frac{\\partial v_y}{\\partial t} = \\frac{\\partial \\sigma_{yx}}{\\partial x} + \\frac{\\partial \\sigma_{yz}}{\\partial z} + f_y. \\notag\n",
    "\\end{equation}\n",
    "\n",
    "We split the estimation of the correct FD operators into 2 steps:\n",
    "\n",
    "1. Where are the material parameters and wavefields on the LHS of the equation located on the staggered grid? In case of the momentum equation they are placed at the gridpoint (i,j).\n",
    "<img src=\"images/SG_SH-Cart.png\" width=\"50%\">\n",
    "2. All derivatives on the RHS have to be centered on the grid point estimated in step 1. Therefore, the spatial derivatives on the RHS of the momentum equation are approximated by  2nd order centered finite-difference operators:\n",
    "\n",
    "\\begin{align}\n",
    "\\biggl(\\frac{\\partial \\sigma_{yx}}{\\partial x}\\biggr)^c(i,j) &\\approx \\frac{\\sigma_{yx}(i+1/2,j) - \\sigma_{yx}(i-1/2,j)}{dx},\\notag\\\\ \n",
    "\\biggl(\\frac{\\partial \\sigma_{yz}}{\\partial z}\\biggr)^c(i,j) &\\approx \\frac{\\sigma_{yz}(i,j+1/2) - \\sigma_{yz}(i,j-1/2)}{dz}.\\notag\n",
    "\\end{align}\n",
    "\n",
    "For the stress-strain relation of the shear stress $\\sigma_{yx}$\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial \\sigma_{yx}}{\\partial t} = \\mu \\frac{\\partial v_{y}}{\\partial x}, \\notag\n",
    "\\end{equation}\n",
    "\n",
    "we use the same approach:\n",
    "\n",
    "1. $\\sigma_{yx}$ on the LHS is located at the grid point (i+1/2,j)\n",
    "<img src=\"images/SG_SH-Cart.png\" width=\"50%\">\n",
    "2. Therefore, we approximate the spatial derivative on the RHS by the 2nd order centered operator:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mu_x \\biggl(\\frac{\\partial v_{y}}{\\partial x}\\biggr)^c(i+1/2,j) \\approx \\mu_x \\frac{v_y(i+1,j) - v_{y}(i,j)}{dx},\\notag \n",
    "\\end{equation}\n",
    "\n",
    "For the stress-strain relation of the shear stress $\\sigma_{yz}$\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial \\sigma_{yz}}{\\partial t} = \\mu \\frac{\\partial v_{y}}{\\partial z}, \\notag\n",
    "\\end{equation}\n",
    "\n",
    "we get:\n",
    "\n",
    "1. $\\sigma_{yz}$ on the LHS is located at the grid point (i,j+1/2)\n",
    "<img src=\"images/SG_SH-Cart.png\" width=\"50%\">\n",
    "2. Therefore, we approximate the spatial derivative on the RHS by the 2nd order centered operator:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mu_z \\biggl(\\frac{\\partial v_{y}}{\\partial x}\\biggr)^c(i,j+1/2) \\approx \\mu_z \\frac{v_y(i,j+1) - v_{y}(i,j)}{dz}.\\notag \n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "So in summary, we get the following spatial FD approximations of the 2D isotropic SH problem:\n",
    "\n",
    "\\begin{align}\n",
    "\\rho \\frac{\\partial v_y}{\\partial t} (i,j)&= \\biggl(\\frac{\\partial \\sigma_{yx}}{\\partial x}\\biggr)^c(i,j) + \\biggl(\\frac{\\partial \\sigma_{yz}}{\\partial z}\\biggl)^c(i,j) + f_y(i,j), \\notag\\\\\n",
    "\\frac{\\partial \\sigma_{yx}}{\\partial t} (i+1/2,j) &= \\mu_x \\biggl(\\frac{\\partial v_{y}}{\\partial x}\\biggr)^c(i+1/2,j),\\notag\\\\\n",
    "\\frac{\\partial \\sigma_{yz}}{\\partial t}(i,j+1/2) &= \\mu_z \\biggl(\\frac{\\partial v_{y}}{\\partial z}\\biggr)^c(i,j+1/2). \\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "Finally, we approximate the time derivatives on the LHS of the partial differential equations using the [**Leapfrog method**](https://en.wikipedia.org/wiki/Leapfrog_integration) which is based on a staggered temporal grid at time steps n and n+1/2:\n",
    "\n",
    "\\begin{align}\n",
    "\\scriptsize\\rho \\frac{v_y(i,j,n+1/2)-v_y(i,j,n-1/2)}{dt}&\\scriptsize= \\biggl(\\frac{\\partial \\sigma_{yx}}{\\partial x}\\biggr)^c(i,j,n) + \\biggl(\\frac{\\partial \\sigma_{yz}}{\\partial z}\\biggl)^c(i,j,n) + f_y(i,j,n), \\notag\\\\\n",
    "\\scriptsize\\frac{\\sigma_{yx}(i+1/2,j,n+1)-\\sigma_{yx}(i+1/2,j,n)}{dt} &\\scriptsize= \\mu_x \\biggl(\\frac{\\partial v_{y}}{\\partial x}\\biggr)^c(i+1/2,j,n+1/2),\\notag\\\\\n",
    "\\scriptsize\\frac{\\sigma_{yz}(i,j+1/2,n+1)-\\sigma_{yz}(i,j+1/2,n)}{dt} &\\scriptsize= \\mu_z \\biggl(\\frac{\\partial v_{y}}{\\partial z}\\biggr)^c(i,j+1/2,n+1/2). \\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "Rearranging to the time steps n+1/2 and n+1 on the LHS, yields the final 2D SH FD scheme on the staggered grid:\n",
    "\n",
    "\\begin{align}\n",
    "\\scriptsize \n",
    "v_y(i,j,n+1/2) &\\scriptsize= \\frac{dt}{\\rho}\\biggl\\{\\biggl(\\frac{\\partial \\sigma_{yx}}{\\partial x}\\biggr)^c(i,j,n) + \\biggl(\\frac{\\partial \\sigma_{yz}}{\\partial z}\\biggl)^c(i,j,n) + f_y(i,j,n)\\biggr\\} + v_y(i,j,n-1/2), \\notag\\\\\n",
    "\\scriptsize\\sigma_{yx}(i+1/2,j,n+1) &\\scriptsize= dt\\; \\mu_x \\biggl(\\frac{\\partial v_{y}}{\\partial x}\\biggr)^c(i+1/2,j,n+1/2) + \\sigma_{yx}(i+1/2,j,n),\\notag\\\\\n",
    "\\scriptsize\\sigma_{yz}(i,j+1/2,n+1) &\\scriptsize= dt\\; \\mu_z \\biggl(\\frac{\\partial v_{y}}{\\partial z}\\biggr)^c(i,j+1/2,n+1/2) + \\sigma_{yz}(i,j+1/2,n), \\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "which we will implement during the next class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "- How to solve the 2D SH problem by the staggered grid finite difference approach"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
